↜ Back to index Introduction to Numerical Analysis 2

Lecture a4

Linear system with a triangular matrix

Given a natural number n \in \mathbb{N} and real numbers a \neq 0 and c, we consider a system of n linear equations of n variables of the form \tag{1} \begin{aligned} a x_1 &= b_1,\\ c x_1 + a x_2 &= b_2,\\ c x_1 + c x_2 + a x_3 &= b_3,\\ &\vdots\\ c x_1 + c x_2 + c x_3 + \cdots + c x_{n-1} + a x_n &= b_n.\\ \end{aligned} For the given right-hand side b_1, b_2, …, b_n, we want to find the values of the unknowns x_1, x_2, …, x_n that satisfy all n equations. We call such x_1, x_2, …, x_n the solution of the linear system.

We can write the same system in a matrix form: \begin{pmatrix} a & 0 &&\cdots & 0\\ c & a & 0 && \vdots\\ c & c & a & 0\\ \vdots& & & \ddots &0\\ c & c & \cdots & c & a\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ \vdots\\ x_n \end{pmatrix} = \begin{pmatrix} b_1\\ b_2\\ b_3\\ \vdots\\ b_n \end{pmatrix} The first matrix is an n\times n square matrix.

The square matrix of the form in the above equation is called lower-triangular matrix (下半三角行列) because all entries above the diagonal are 0.

Example 1. Let n = 3, a = 2 and c = 3. Then the system reads as \begin{aligned} 2 x_1 &= b_1,\\ 3 x_1 + 2 x_2 &= b_2,\\ 3 x_1 + 3 x_2 + 2 x_3 &= b_3.\\ \end{aligned} and in matrix form \begin{pmatrix} 2 & 0 &0\\ 3 & 2 & 0\\ 3 & 3 & 2 \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix} = \begin{pmatrix} b_1\\ b_2\\ b_3 \end{pmatrix} We can find the solution easily as \begin{aligned} x_1 &= \tfrac 12 b_1,\\ x_2 &= \tfrac 12 (b_2 - 3 x_1),\\ x_3 &= \tfrac 12 (b_3 - 3 x_1 - 3 x_2).\\ \end{aligned}

Here is an example of a Fortran code that computes the solution:

implicit none
integer, parameter :: n = 3
real b(n), x(n)

print *, "Enter b(1), b(2) and b(3)"
read *, b(1), b(2), b(3)

x(1) = 0.5 * b(1)
x(2) = 0.5 * (b(2) - 3 * x(1))
x(3) = 0.5 * (b(3) - 3 * x(1) - 3 * x(2))

print *, "solution = "

print *, x(1), x(2), x(3)

end

We can test it easily:

$ gfortran example.f90
$ ./a.exe
 Enter b(1), b(2) and b(3)
1 2 3
 solution = 
  0.500000000      0.250000000      0.375000000

Replace ./a.exe with ./a.out if you are on macOS or Linux.

Testing tip. To not have to enter the numbers manually every time we run the program, we can use the following trick using the printf command and the pipe operator |:

$ printf "1 2 3" | ./a.exe
 Enter b(1), b(2) and b(3)
 solution = 
  0.500000000      0.250000000      0.375000000

or even shorter

$ ./a.exe <<< "1 2 3"
 Enter b(1), b(2) and b(3)
 solution = 
  0.500000000      0.250000000      0.375000000

This way it is easy to re-run the code with the same input many times.

Finally, it is also possible to read the input data from a file instead. First, in Emacs create file input.txt (the name is not important) with the content:

1
2
3

(On the command line you can use printf "1\n2\n3" > input.txt.)

Then run the code as:

$ ./a.exe < input.txt
 Enter b(1), b(2) and b(3)
 solution = 
  0.500000000      0.250000000      0.375000000

Report 1 Problems

Please follow the instructions:

  1. For each of the following problems, submit the original fortran source code. Do not submit Emacs’ backup files like name.f~ or name.f90~ ending with ~!
  2. Make sure that the program compiles using gfortran.
  3. After it compiles, make sure that the output of the program matches the output given in the attached example when run the way the example indicates. In particular, there should be no extra things in the output like intermediate results.
  4. Upload the code to Acanthus portal at the appropriate problem at: by 2021/11/5 (Fri).

You can get up to 20 points for each problem. The full score is 100 points.

Important. Upload only your own code; do not show your code to others. I check for similar code.

Problem 1.

Write a Fortran program that reads a and b from the user (standard input) and prints the solutions x of the linear equation a x = b.

Make sure to check that a \neq 0.

a=2, b = 3 \Rightarrow x = \tfrac 32

$ printf "2\n3" | ./a.exe
   1.50000000

If a = 0 the program should print a message.

$ printf "0\n3" | ./a.exe
 Error: a must be nonzero

Problem 2.

Write a Fortran program that reads a, c and b_1, b_2, …, b_5 from the user (from the standard input) and prints out the solution x_1, x_2, …, x_5 of the linear system (1) with n = 5. The program may assume that a \neq 0 (no need to check).

Hint. I recommend that you use a do loop.

a=2, c = 5, b = (-1, 2, 2, 3, 2)

$ printf "2\n5\n-1\n2\n2\n3\n2" | ./a.exe 
 -0.500000000       2.25000000      -3.37500000       5.56250000      -8.84375000

Problem 3. Modify the program in Problem 2 so that it reads n followed by a, c and b_1, b_2, …, b_n and prints out the solution x_1, x_2, …, x_n of the linear system (1).

n=3, a = 2, c = 3, b = (1, 2, 3) \Rightarrow x = (\tfrac 12, \tfrac 14, \tfrac 38)

$ printf "3\n2\n3\n1\n2\n3" | ./a.exe 
  0.500000000      0.250000000      0.375000000

n = 10, a = 1, c = 2, b = (1, 2, \ldots, 10)

$ printf "10\n1\n2\n1\n2\n3\n4\n5\n6\n7\n8\n9\n10" | ./a.exe
   1.00000000       0.00000000       1.00000000       0.00000000       1.00000000       0.00000000       1.00000000       0.00000000       1.00000000       0.00000000

Problem 4. Modify the program in Problem 3 so that it prints out the solution of the system with the upper triangular matrix instead: \begin{pmatrix} a & c &c&\cdots & c\\ 0 & a & c && c\\ & 0 & a & &\vdots\\ \vdots& & & \ddots &c\\ 0 & \cdots & 0 & 0 & a\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ \vdots\\ x_n \end{pmatrix} = \begin{pmatrix} b_1\\ b_2\\ b_3\\ \vdots\\ b_n \end{pmatrix}

n=3, a = 2, c = 3, b = (1, 2, 3) \Rightarrow x = (\tfrac 18, -\tfrac 54, \tfrac 32)

$ printf "3\n2\n3\n1\n2\n3"  | ./a.exe
  0.125000000      -1.25000000       1.50000000

An example with n=5

$ printf "5\n2\n1\n2\n2\n3\n4\n1" | ./a.exe
 -0.156250000     -0.312500000      0.375000000       1.75000000      0.500000000

Problem 5. Combine the programs in Problem 3 and 4 so that the resulting program reads n , c, d and b_1, b_2, …, b_n and prints the solution x_1, x_2, …, x_n of the linear system: \begin{pmatrix} 1 & 0 &&\cdots & 0\\ c & 1 & 0 && \vdots\\ c & c & 1 & 0\\ \vdots& & & \ddots &0\\ c & c & \cdots & c & 1\\ \end{pmatrix} \begin{pmatrix} 1 & d &d&\cdots & d\\ 0 & 1 & d && d\\ & 0 & 1 & &\vdots\\ \vdots& & & \ddots &d\\ 0 & \cdots & 0 & 0 & 1\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ \vdots\\ x_n \end{pmatrix} = \begin{pmatrix} b_1\\ b_2\\ b_3\\ \vdots\\ b_n \end{pmatrix}

Hint. Note that MN x = b is equivalent to M(Nx) = b and so we first solve My = b for y and then Nx = y for x.

n=3, c = 1, d = 2, b = (1, 2, 3) \Rightarrow x = (1, -1, 1)

$ printf "3\n1\n2\n1\n2\n3" | ./a.exe
   1.00000000      -1.00000000       1.00000000

n=8, c = 2, d = 3, b = (1, 2, 3, 2, 1, 2, 3, 4)

$ printf "8\n2\n3\n1\n2\n3\n2\n1\n2\n3\n4" | ./a.exe
   151.000000      -75.0000000       37.0000000      -17.0000000       7.00000000      -3.00000000       1.00000000       0.00000000